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Correlator product states (CPS) are a class of tensor network wavefunctions applicable to strongly 
correlated problems in arbitrary dimensions. Here, we present a method for optimizing and eval- 
uating the energy of the CPS wavefunction that is non-variational but entirely deterministic. The 
fundamental assumption underlying our technique is that the CPS wavefunction is an exact eigen- 
state of the Hamiltonian, allowing the energy to be obtained approximately through a projection of 
the Schrodinger equation. The validity of this approximation is tested on two dimensional lattices 
for the spin-i antiferromagnetic Heisenberg model, the spinless Hubbard model, and the full Hub- 
bard model. In each of these models, the projected method reproduces the variational CPS energy 
to within 1%. For fermionic systems, we also demonstrate the incorporation of a Slater determi- 
nant reference into the ansatz, which allows CPS to act as a generalization of the Jastrow-Slater 
wavefunction. 



I. INTRODUCTION 

Tensor network wavefunctions represent an efficient 
approach to modeling strongly interacting quantum sys- 
tems. While the density matrix renormalization group 
(DMRG) algorithm [US] and corresponding matrix prod- 
uct state (MPS) have been extremely successful in one 
dimensional systems, the application of tensor networks 
to two dimensions remains a work in progress. The 
pair entangled product state (PEPS) wavefunction 0- 
|l3j represents the natural extension of the MPS in two 
dimensions, but lacks the tractability with respect to 
evaluating expectation values that the MPS enjoys in 
one dimension. Similarly, the multi-scale entanglement 
renormalization ansatz (MERA) [l4], [H| has been lim- 
ited by computational cost. Recently, a simpler class 
of tensor network wavefunction has been proposed by 
both some of the present authors [l6| and others [171 - 
[l9l | . This wavefunction, known as the correlator product 
state (CPS), the entangled plaquette state (EPS), or the 
complete graph tensor network, seeks to correlate neigh- 
boring lattice sites directly, without the use of auxiliary 
bond indices. It is closely related to the string bond 
state, [13, [2l[ although the correlators are simple ten- 
sors and do not involve any underlying MPS structure. 
As discussed previously [lfj , the CPS ansatz also encom- 
passes a number of other wavefunctions, including the 
Huse-Elser ansatz I22l. the Laughlin wavefunction [23| . 
and the toric code [241 ] . So far in its development, the 
CPS ansatz has been studied through the use of vari- 
ational Monte Carlo (VMC) techniques, which provide 
efficient procedures for its optimization and the evalua- 
tion of expectation values in two and higher dimensions 

Here we present an alternative approach to optimizing 
and evaluating the energy of the CPS wavefunction that 
is entirely deterministic, i.e. no stochastic sampling is 
involved. Much like the coupled cluster (CC) formalism 
commonly used in quantum chemistry [251 ] . our method 



relies on assuming that the CPS wavefunction is an ex- 
act eigenstate of the Hamiltonian. Under this assump- 
tion, the Schrodinger equation may be projected against 
a carefully chosen bra state in order to yield an approxi- 
mate, non-variational expression for the system's energy 
that can be evaluated without stochastic sampling. To 
be accurate, this projected energy functional requires an 
optimization of the CPS wavefunction that differs from 
that used in VMC. As the accuracy of the functional de- 
pends on how close the CPS state is to a true eigenstate, 
the wavefunction is optimized by satisfying a set of pro- 
jected Schrodinger equations, which comprise a subset 
of the conditions necessary for a wavefunction to be a 
Hamiltonian eigenstate. Central to our approach is the 
formulation of the CPS wavefunction as a product of lo- 
cal, invertible operators acting on a reference wavefunc- 
tion. In previous work, the CPS reference has been taken 
as an equally weighted sum of the states in the Hilbert 
space, which we name the uniform reference. Here, we 
also consider a Slater determinant reference for fermionic 
systems, which allows CPS to act as a generalization of 
the Jastrow-Slater wavefunction long used in variational 
Monte Carlo simulations [26l.[27j. 

To evaluate the accuracy of our projected energy func- 
tional, we have applied both the VMC CPS and the pro- 
jected CPS (PCPS) methods to three model Hamiltoni- 
ans on two-dimensional lattices. Our first test case is the 
spin-i antiferromagnetic Heisenberg model, which is cur- 
rently the most widely studied model for the CPS ansatz. 
[pol [ItJ We then tested our method on two fermionic 
models, the spinless and regular Hubbard models, in 
which we employed the Slater determinant reference. In 
all cases, we find that the PCPS method reproduces the 
results of VMC CPS to within 1% in the ground state 
energy. This represents a first proof of principle that the 
technique of projecting the Schrodinger equation can be 
effectively applied to a tensor network wavefunction. 

We begin by reviewing the structure of the CPS ansatz 
and presenting its formulation in terms of operators act- 
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ing on a reference wavefunction (Sec. Ill A|) . We then 
review how the energy may be evaluated stochastically 
through VMC fSec. lIIB] ) before presenting the projected 
energy functional (Sec. Ill C[) and how a system of pro- 
jected Schrodinger equations may be used to optimize 
the CPS ansatz (Sec. Ill D|) . Results are then presented 
for the Heisenberg model (Sec. IIII A"|) . spinless Hubbard 
model (Sec. HUH), and full Hubbard model fSec. iHTCl) . 
Finally, we summarize the main points and discuss future 
directions for CPS wavefunction research fSec. IIV|) . 



II. THEORY 



The CPS Ansatz 



Central to the definition of the CPS wavefunction is 
the concept of a correlator, which directly encodes cor- 
relations between lattice sites. Consider a subset of the 
lattice sites, P = {i\,i-2 ■ ■ -ii}, where each site i carries 
quantum states |rij) such as spin or fermionic degrees 
of freedom. A correlator with domain P is defined by 
its amplitudes Cp 1 ' 2 " H for each configuration in the 
Hilbcrt space of P. The total CPS wavefunction \^>) is 
built from multiple correlators on separate (possibly over- 
lapping) sets of the sites. In a simple CPS, for a lattice 
configuration \n\n 2 ■ ■ ■ n k ), the corresponding wavefunc- 
tion amplitude (nin2 ■ ■ ■ n k \^) is given by the product of 
correlator amplitudes, i.e. 

I*>= E Y[Cp in ^- nH \n in2 ...n k ). (1) 
n1n2--.ru P 

For example, for nearest neighbor correlators in one- 
dimension with open boundary conditions, Eq. ([!]) is 
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z \nin 2 ...n k ). (2) 



The CPS wavefunction can also be written in a more 
general way so that the correlators appear as operators 
acting on a reference state. We first define the projection 
operator p rii = \rii)(m\ which projects to the subspace in 
which site i is in state n$. For example, 



p„> |ni . . 



■ ■ n k ) = 8n] \ni...m... n k ), 



(3) 



where 5 is the Kronecker delta. The correlator opera- 
tor for domain P — i 2l i{\ can now be defined by 
combining the correlator amplitudes with a product of 
projection operators for each site in P, 
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(4) 



We then write the CPS wavefunction as a product of 
correlator operators acting on a reference state |$), 



|¥) = IJ6 P |$). 



(5) 



In order for this wavefunction to be equivalent to that in 
Eq. ([1]), we must choose |<&) to be the sum of all possible 
lattice configurations with equal coefficients, a state we 
refer to as the uniform reference, 
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} J \nin 2 .-.n k ). 
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n 1 n 2 ...n k 



In this work, we usually restrict the summation in Eq. 
^ to states with particular quantum numbers. For ex- 
ample, for a spin system, we sum only over those spin 
configurations with a given S z value, while for a fermionic 
system, we sum only over configurations with a given N 
and given S z . 

One can also use other reference functions for |$), in 
which case the correlator operators can be seen as provid- 
ing some additional correlation beyond what is present 
in the reference. In practice, only some reference func- 
tions will be useful. Specifically, they must not interfere 
with the efficient evaluation of expectation values. In 
this work, we will also explore the Slater determinant as 
a reference function, which allows the CPS wavefunction 
to act as a generalization of the Jastrow- Slater wavefunc- 
tion [IE 51. 



B. The Monte Carlo Variational Energy 



In previous work, |16l - ll8| the CPS energy has been 
evaluated by means of Monte Carlo (MC) sampling. For 
a general wavefunction, 

l*)= E * nin2 - nk \nin 2 ...n k ), (7) 

n 1 n 2 ...n k 

the energy may be written as 

<?s\Hm 
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where the local energy Ei J (nin 2 ...n k ) is defined by 



E h (mn 2 ...n k ) 



tyn 1 n 2 ...n h 



(n 1 n 2 ...n k \H\n[n' 2 ...n' k ). (9) 



A Markov chain can then be used to sample the prob- 
ability distribution \^nm 2 ...n k |2/(^|^ anc i compute the 
overall energy as the average of the sampled local en- 
ergies. As discussed in previous reports 0113,1131, the 
CPS wavefunction is well suited to MC sampling because 



the amplitudes \$ n i n 2--- n k are eaS y to evaluate. When |$) 
is chosen as the uniform reference, the amplitudes are 
simply a product of correlator values, as in Eq. (p}. When 
|$) is allowed to take on other forms, the contribution of 
the reference function, the factor (niri2...nk\'&), must be 
included as well. In order to evaluate the energy by MC 
sampling in practice, we are therefore limited to forms of 
|$) for which this factor can be evaluated efficiently. Two 
forms of |$) have long been used in quantum Monte Carlo 
simulations. For the Slater determinant, (niri2---nk\&) is 
equal to a determinant of orbital coefficients that can be 
evaluated in at most O(Np) time by an LU decomposi- 
tion, where N p is the number of particles in the system. 
In practice, this cost can be reduced to 0{N p ) during a 
Markov chain iteration by using the matrix determinant 
lemma and the Sherman-Morrison formula [28l 29(. An- 
other reference function for which sampling is efficient 
is the BCS wavefunction [30| and its number projected 
form, the antisymmetrized geminal power [3 lM33j ] . How- 
ever, we will not explore the BCS reference here. 



C. The Deterministic Projected Energy Functional 

In the current work, in addition to using a Monte Carlo 
evaluation of the energy, we also explore another method 
which does not require the use of stochastic sampling. We 
first make an approximation in which we assume that the 
CPS wavefunction is an eigenstate of the Hamiltonian, 



o o o o o o 
o o o o o o 

O O © <8> O O 

o o o o o o 
o o o o o o 

FIG. 1: An example 2D lattice showing the cluster (shaded) 
resulting from nearest neighbor pair correlators around a hop- 
ping operator. The site touched by the creation operator is 
marked with a "+" and the site touched by the destruction 
operator by an "x". Computing the projected CPS energy 
requires summing over only the configurations of the differ- 
ent clusters produced by the Hamiltonian's operators, rather 
than over all lattice configurations. The energy for a local 
Hamiltonian can thus be evaluated in 0(Nd M ) time, where 
N is the size of the lattice, d the number of configurations of 
a site, and M the number of sites in a cluster. 



(10) 



We then define an inverse correlator product bra state 
(^|, which is obtained by left multiplying (<I>| with a 
product of inverse correlator operators 



where 



i 



(ii) 



„ni 1 «i 2 ...n-i Prii 1 Prii 2 ■ • ■ * (12) 
ni 1 rii 2 ...rii l r 



This gives an approximate energy functional that is exact 
if |^>) is an eigenstate, 



This hopping operator is associated with two sites x and 
y. We can divide correlators (and their inverse correla- 
tors) into two kinds: those which touch the sites x or y 
and those which do not. Correlators which do not touch 
x or y can be commuted past the the hopping operator 
and thus cancel with their inverse partners. The set of 
correlators which touch x or y define a connected cluster 
of sites, fl xy = {ci, C2, . . . , c xy } (see Fig. [T] for an exam- 
ple). The central point is that, as long as the domain size 
of the correlators is independent of lattice size, the size 
of the cluster Q xy is also independent of lattice size, and 
this allows expectation values of the form (|T4| to be eval- 
uated efficiently (i.e. with a cost scaling polynomially in 
the lattice size). 



e = mmm) 



(13) 



Provided that the correlators are not too large or long- 
ranged, Eq. (jT5)) can be evaluated efficiently by consid- 
ering the different operators inside the Hamiltonian in- 
dividually. Consider the energy contribution E xy from a 
hopping operator a x a y , 

s xy ~ (*|4<g*) = (i'inc Q 1 «kn^l*>- (i4) 

Q P 



Commuting terms past a x a y , the product of correlators 
in (fl4|) becomes 



E xy = m II Cq 1 ala y J] Cp\$), (15) 



where S xy denotes the set of correlators that touch x or y. 
For a more explicit expression, we separate the correla- 
tors into their amplitude and projection operator compo- 
nents. The expectation value of the projection operators 
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define a many-body reduced density-matrix (RDM) 7, 
%,° 1 '.'.2' ! ° a = ■ • ■ Pnjn y ■ ■ ■ /3„ 

4««Pn' cl ■ ■ • Pn'Jn' y ■ ■ ■ pn' Cxy |$)- (16) 

Combining this density matrix with the correlator am- 
plitudes, we obtain 

E ** = E ^;::. n > n §W ( 17 ) 

n' ...n' 

c l c :ry 

For each term in the sum, the quotient of correlators is 
easily evaluated. However, we must also evaluate the rel- 
evant RDM element, which limits us to references for 
which 7 is readily available. The uniform reference's 
RDM elements are trivial to evaluate (even with par- 
ticle number and S z restrictions) while those of a Slater 
determinant can be computed as a determinant of one- 
body RDM elements in N 3 time, where N is the number 
of lattice sites. With the RDM elements in hand, the 
final task is to sum over all terms in Eq. (fTT|) . Due to 
the sparsity of 7, there are at most d M non-zero terms, 
where d is the size of the single-site Hilbert space and M 
is the number of sites in the cluster. For sufficiently small 
and local correlators (i.e. for small enough clusters), the 
summation can be carried out exactly without resorting 
to importance sampling. To summarize, each operator 
within the Hamiltonian defines a cluster, on which we 
compute a sum of correlator ratios and RDM elements. 
The total energy is then the sum of each operators' con- 
tribution. 

Although the energy functional of Eq. (JT3J) allows for a 
non-stochastic evaluation of the CPS wavefunction's en- 
ergy, it possesses a number of shortcomings that should 
be discussed. First, it is not a variational form, as the ini- 
tial assumption that the wavefunction is an exact Hamil- 
tonian eigenfunction is only an approximation. Second, 
the functional is only tractable for small, local correla- 
tors. For example, if all lattice site pairs were used to de- 
fine correlators, then each hopping operator would create 
a cluster the size of the lattice, and the functional would 
be no more tractable than the variational energy expres- 
sion. Similarly, local correlators that contain a large 
number of sites may lead to clusters whose configurations 
are too numerous to sum over explicitly. A final concern 
with the energy functional is the form of the Hamilto- 
nian. In lattice models, one rarely encounters nonlocal 
operators or operators that touch more than two sites, 
so the clusters formed around them will be manageable 
provided the correlators are small and local. However, in 
quantum chemistry, the Hamiltonian contains the long 
range Coulomb interaction, which manifests as a nonlo- 
cal four-site operator. This would severely restrict the 
size of correlator for which summation over the cluster 
would be feasible. 



D. Optimizing the Wavefunction by Projection 

When using the variational Monte Carlo energy func- 
tional, we can obtain the ground state by minimizing 
the energy. Several different algorithms can be used. In 
Ref. we used a "sweep" algorithm where a single 

correlator is updated at a time by solving an effective 
Schrodinger equation. For the variational Monte Carlo 
results reported in this work, we have employed a steepest 
descent optimization scheme that updates all correlators 
simultaneously. 

The deterministic projected energy functional, how- 
ever, is not a variational expression. Consequently we 
require a different method to optimize the correlators 
in the CPS when using this functional. Indeed, accu- 
rate results from the projected functional require corre- 
lators that make the wavefunction as close as possible to 
a Hamiltonian eigenstate. This optimization condition is 
slightly different from minimizing the energy, and will in 
general produce slightly different correlators. 

We begin as in the previous section by assuming that 
the CPS wavefunction is a true Hamiltonian eigenstate 
and therefore satisfies the Schrodinger equation, 

(H-E) J[Cp\<$>} ~0. (18) 
p 

We then seek to enforce this condition approximately by 
applying a number of projections. We define a set of bra 
states (fyp 1 12 H \ obtained from the inverse CPS bra 
state in Eq. (fTTj) 

{^ P iini2 '" ni, \ = {^\p nil p ni2 ...p nil : h,i 2 -.-ii e P, 

(19) 

where the projectors are those associated with a given 
correlator Cp. Projecting with these bra states, we ob- 
tain 

A n« 1 n ta ...n t| = ^n^...n H ^ _ £ ^ _ Q (20) 

By requiring that each ^i"' 2 "'"'i vanish, we create a 
system of nonlinear projected Schrodinger equations that 
can be solved for the correlator amplitudes. As with the 
energy, we can evaluate the contributions of the individ- 
ual Hamiltonian elements to these equations separately. 
The local cluster formed around each Hamiltonian op- 
erator is the same as in the energy evaluation, except 
that the projection operators p ni fix some of the sites' 
occupations. 

In order to solve the system of equations given in Eq. 
(l20l) . we have taken two approaches. First, the deriva- 
tives of each A p n * 2 " H with respect to each correla- 
tor amplitude can be evaluated, producing a Jacobian 
matrix that can be used in a standard Newton-Raphson 
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TABLE I: Energies per site of the antiferromagnetic spin-^ Heisenberg model with periodic boundary conditions and local 
correlators. The correlator sizes are given with the method names. Energies are in units of J, with the number in parentheses 
representing the uncertainty in the final digit. See Sec. IIII Al for details. 



Lattice Size 


PCPS 2x2 


VMC CPS 2x2 


VMC CPS 3x3 


VMC EPS 4x4 a 


SSE b 


4x4 


-0.694226 


-0.69432(1) 


-0.70150(1) 


-0.7016(1) 


-0.701777(7) 


6x6 


-0.668615 


-0.66806(2) 


-0.67615(1) 


-0.6785(2) 


-0.678873(4) 


8x8 


-0.665920 


-0.66165(1) 


-0.66989(1) 


-0.6724(3) 


-0.673487(4) 


10x10 


-0.664910 


-0.65900(1) 


-0.66776(1) 


-0.6699(3) 


-0.671549(4) 



a Entangled plaquette state, see Ref. [171 . 
b Stochastic series expansion, see Ref. [341 . 



procedure. Alternatively, we may take the approach of 
constructing and diagonalizing a local Hamiltonian for 
each correlator. This second approach arises from fixing 
all variables except for one correlator's amplitudes, which 
makes the equations linear in those amplitudes. The ex- 
act solution to this simplified system of linear equations 
is then equivalent to the lowest eigenvector of the local 
Hamiltonian. The local Hamiltonian for a correlator on 
sites P = i2--.ii} is defined by 

-i> II <V /'<> •••/'" (2i) 

1 *' Q^p 

{H-E)p n ^...p Ki l[ C Q ,\<t>). 

As before, we may evaluate this expression efficiently 
if we consider each Hamiltonian operator's contribution 
separately, in which case we need only sum over the 
configurations of the local cluster around each opera- 
tor. Once we have evaluated the local Hamiltonian and 
found its lowest eigenvalue and eigenvector, we update 
the correlator corresponding to P with the eigenvector's 
elements. This process is performed for each correla- 
tor, and then is repeated iteratively until A is sufficiently 
small. 



III. RESULTS 



the projected method (see Sec. Ill C|) restricted us to us- 
ing 4-site (2x2) correlators. The variational Monte Carlo 
method can practically support up to 16 site (4x4) cor- 
relators, results for which were reported earlier by Mez- 
zacapo et al [lj|- We have reproduced the 4x4 results 
and also optimized the variational wavefunction for 4 site 
(2x2) and 9 site (3x3) correlators. The results are sum- 
marized in Table U for L = 4, 6, 8, and 10. For 2x2 
correlators, both the variational and projected methods 
produce energies with a relative error between 1 and 2% 
of the essentially exact stochastic series expansion (SSE) 
[3~il ]. Thus we see that the eigenfunction assumption un- 
derpinning the projected method is already reasonable 
for 2x2 correlators. We would naturally expect the pro- 
jected and VMC methods to produce even closer energies 
to each other for larger correlators, as in this case the 
wavefunction would be even closer to an exact Hamilto- 
nian eigenstate. By using larger correlators in the vari- 
ational wavefunction, we see that the relative error can 
be systematically reduced. The 3x3 correlators produce 
relative errors below 0.6%, while the 4x4 correlators pro- 
duce relative errors below 0.25%. In summary, the pro- 
jected method produces reasonably accurate results on 
the 2D Heisenberg model that are comparable to those 
of the variational method, however the limitation on cor- 
relator size prevents it from obtaining the higher accu- 
racies achieved with large correlators in the variational 
method. 



A. Antiferromagnetic Heisenberg Model 

We have optimized the CPS ansatz using the vari- 
ational Monte Carlo and projected energy functionals, 
for the antiferromagnetic spin-^ Heisenberg model on a 
number of periodic square lattices of edge length L. The 
Hamiltonian, in which the parameter J is assumed to be 
positive, is written as 

H = S-j, (22) 

where the notation indicates the set of all nearest 
neighbor pairs. The limitations of the cluster size for 



B. Spinless Hubbard Model 

We have treated a 20-site (4x5) spinless Hubbard lat- 
tice with open boundary conditions using 4 site (2x2) 
correlators and the Hartree-Fock Slater determinant ref- 
erence. While the CPS ansatz is capable of treating much 
larger lattices, we chose this lattice size in order to com- 
pare to exact results. The Hamiltonian for the spinless 
Hubbard model is 

H = y~] Ua}aia^aj — t(a}aj + a^ai), (23) 
(i,3) 
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n 1 1 1 1 1 1 r 



Hartree-Fock — 0— 
VMC CPS 2x2 e 
PCPS2x2 — •- 
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FIG. 2: Ground state energy errors for the 20 site (4x5) 
spinless Hubbard model with 10 particles and open bound- 
ary conditions. Exact results were computed using the ALPS 
program [35l ]. 



l 1 1 1 1 1 1 1 r 



Hartree-Fock — 0-- 

VMC CPS 2x2 a- 

PCPS2x2 — •- 




FIG. 3: Ground state energy errors for the 20 site (4x5) spin- 
less Hubbard model with 9 particles and open boundary con- 
ditions. Exact results were computed using the ALPS pro- 
gram [35l ]. 



in which a] and are the fermionic particle creation and 
destruction operators on site i, and represents the 
sum over all nearest neighbor site pairs. The first term in 
Eq. (l23|) represents a repulsive interaction between par- 
ticles on neighboring sites, while the second term repre- 
sents hopping between sites. 

The quality of the projected energy functional is sensi- 
tive to the quality of the reference function |<£>). Using the 
uniform reference is therefore undesirable, as it is clearly 
incorrect for large U/t. We instead use a Slater deter- 
minant whose orbital coefficients have been variationally 
optimized through Hartree-Fock theory. Such a determi- 
nant satisfies the model exactly in the small U/t limit 
for any filling, and in the large U/t limit at half filling. 
Furthermore, we have observed that the determinant is 
also accurate for single hole doping at large U /t. 

Results for half-filling and single hole doping are pre- 
sented in Figs. [2] and El respectively. For the case of 
half-filling, the CPS results using both the projected and 
variational functions are within 1% of the exact energy 
for all values of the ratio U /t. For single hole doping, the 
worst error is below 3%. If one considers the CPS wave- 
function as a generalization of the Jastrow- Slater wave- 
function, these accuracies are not extraordinary. How- 
ever, it should be stressed that the projected CPS re- 
sults, which are essentially the same as the variational 
results, require no stochastic sampling as is typically used 
in Jastrow-Slater simulations, but are obtained entirely 
deterministically. 



C. Full Hubbard Model 

We have applied the CPS ansatz to the Hubbard model 
with open boundary conditions and half filling in both 
one and two dimensions. The Hubbard Hamiltonian is 

H = U^2 a W a li a 4 " t ( a h a j° + a i<r a ^)> 

(24) 

in which o|f(,|.) and a,if(i) are the fermionic creation and 
destruction operators for particles with spin t (4-), and 
the notation refers to the set of all nearest neigh- 
bor site pairs. We use as our reference function the re- 
stricted Hartree-Fock (RHF) Slater determinant. While 
this reference is far from ideal as it is qualitatively in- 
correct in the large U/t limit, it is sufficient for illustrat- 
ing the central point we seek to convey. Better energies 
could be obtained with an unrestricted determinant or a 
particle number projected BCS reference [3fl l3lj|. but in 
the present discussion we limit ourselves to asking how 
the projected and variational CPS formulations compare, 
and for this purpose the RHF determinant is sufficient as 
a reference function. 

Results for the ratios U/t = 2 and U/t — 4 are pre- 
sented in Table Ull For the lower ratio, the RHF reference 
produces energies in error by 6-10%, which are reduced to 
1% or less after the optimization of either the projected 
or variational CPS wavefunction. More interesting is the 
fact that the projected and variational CPS energies dif- 
fer from each other by less than 0.02% for both the one 
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TABLE II: Ground state energies for the Hubbard model at half filling with open boundary conditions. The DMRG results 
used m=1600 renormalized states. The correlator sizes used were 3-site for the ID lattices and 4-site (2x2) for the 4x5 lattice. 
Energies are in units of t, with the uncertainty in the final digit placed in parentheses. See Sec. IIII Cl for details. 

U/t = 2 



Lattice Size 



DMRG 



RHF 



PCPS 



VMC CPS 



1x14 
1x18 
1x22 
4x5 



-11.279897 
-14.653987 
-18.029379 
-20.127521 



-10.133544 
-13.219131 
-16.307287 
-18.800678 



-11.241776 
-14.592961 
-17.946260 
-19.920320 



-11.240(1) 
-14.591(1) 
-17.947(2) 
-19.917(1) 



U/t = 4 


Lattice Size 


DMRG 


RHF 


PCPS 


VMC CPS 


1x14 


-7.672349 


-3.133544 


-7.631100 


-7.556(1) 


1x18 


-9.965398 


-4.219131 


-9.842409 


-9.770(3) 


1x22 


-12.259082 


-5.307287 


-12.042344 


-11.964(4) 


4x5 


-14.404488 


-8.800678 


-13.384297 


-13.350(1) 



and two dimensional lattices. For the case of U/t = 4, 
the RHF reference is qualitatively incorrect with rela- 
tive errors as high as 60%. Despite this poor starting 
point, both variational and projected CPS reduce the 
error to 2% and 7% in one and two dimensions, respec- 
tively. These accuracies are worse than in the U/t = 2 
case, but this is primarily due to the poor reference and 
could be improved. Crucially, the projected method still 
effectively reproduces the variational CPS energy, the rel- 
ative difference between them never exceeding 1%. 

IV. CONCLUSIONS 

We have shown that by applying certain projections 
to the Schrodinger equation, the correlator product state 
(CPS) wavefunction can be optimized and its energy ap- 
proximated without the need for stochastic sampling. 
The energies produced by this projected CPS method 
differ by less than 1% from the corresponding variational 
Monte Carlo (VMC) results in tests on three types of 
two dimensional systems: the spin-^ antiferromagnetic 
Heisenberg model, the spinless Hubbard model, and the 
full Hubbard model. While the projection procedure is 
currently limited to smaller correlators than the VMC ap- 



proach, it is nonetheless encouraging that the projected 
Schodinger equation technique can be applied to a tensor 
network wavefunction, and it is our hope that this tech- 
nique will find use with other tensor network ansatze. 

We have also demonstrated that the CPS wavefunc- 
tion may be usefully separated into a product of corre- 
lator operators acting on some reference function. For 
fcrmions, we have showed that reasonable results can be 
obtained through the use of a Slater determinant refer- 
ence, which makes the CPS ansatz a generalization of the 
Jastrow-Slater wavefunction. For the Hubbard model, 
we employed a restricted determinant, and thus we ex- 
pect further improvements in accuracy could be achieved 
by simply breaking the spin restriction. Other candi- 
dates for use as CPS references include the BCS and 
AGP wavefunctions, which should further improve on the 
Slater determinant and will be explored in future work. 
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